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Eigenmode Analysis of Ducted Flows With 
Radially-Dependent Axial and Swirl Velocity Components 


Summary 


This report characterizes the sets of small disturbances possible in cylindrical and annular ducts with mean 
flow whose axial and tangential components vary arbitrarily with radius. The linearized equations of motion 
axe presented and discussed, and then exponential forms for the axial, circumferential, and time dependencies 
of any unsteady disturbances 'are assumed. The resultant equations form a generalized eigenvalue problem, 
the solution of which yields the axial wavenumbers and radiad mode shapes of the unsteady disturbances. 
Two numerical discretizations are applied to the system of equations: (1) a spectral collocation technique 
based on Chebyshev polynomial expansions on the Gauss-Lobatto points, and (2) second and fourth order 
finite differences on uniform grids. The discretized equations are solved using a standard eigensystem package 
employing the QR algorithm. 

The eigenvalues fall into two primary categories: a discrete set (analogous to the acoustic modes found 
in uniform mean flows) and a continuous band (analogous to convected disturbances in uniform mean flows) 
where the phase velocities of the disturbances correspond to the local mean flow velocities. Sample mode 
shapes and eigensystem distributions are presented for both sheared axial and swirling flows. 

The physics of swirling flows is examined with reference to hydrodynamic stability and completeness 
of the eigensystem expansions. The effect of assuming exponential dependence in the axial direction is 
discussed. 
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1. Introduction 


Our understanding of ducted, unsteady turbomachinery flows depends to a considerable extent on un- 
derstanding the behavior of small disturbances in the fluid. For uniform mean flows, these disturbances 
can be classified as vortical and entropic waves, which are purely convected by the mean flow, and acoustic 
waves, which either propagate unattenuated or decay exponentially away from their source of origin [1]. The 
response to any arbitrary unsteady excitation can then be determined as a linear combination of responses 
due to the entropic, vortical, and acoustic parts of the excitation. 

For uniform mean flows, the acoustic disturbances in an annular or cylindrical duct cam be described 
in terms of individual modes [2]. The modes have exponential dependencies in time and in the axial and 
circumferential directions. The radial behavior is then described in terms of Bessel functions of the first 
and second kind. The propagation characteristics of the acoustic modes are dependent upon the mean 
flow quantities, duct geometry, and unsteady excitation parameters. By an appropriate selection of these 
quantities, undesired acoustic modes can be “cut off,” i.e., generated such that they decay exponentially as 
they move in the axial direction. 

The unsteady modes also serve as means of communication between acoustic elements in a turbomachine 
[3, 4], and can be used to determine the coupled responses of several acoustic'elements simultaneously [5]. In 
addition, the derivation of so-called “non-reflecting” far-field boundary conditions for CFD analyses depends 
on an accurate modeling of the far-field unsteady behavior, otherwise the inlet and exit boundaries will give 
rise to spurious, non-physical reflections [1, 6, 7]. 

The behavior of unsteady perturbations in ducted, uniform axial mean flows can be described by a 
convected wave equation. For uniform mean flows, this equation is separable into the dependencies listed 
above. To determine the analogous behavior for nonuniform mean flows, the assumption of periodic behavior 
in time and in the axial and tangential directions results in a linear second-order ordinary differential equation 
for the radial behavior. This equation, while linear, has variable coefficients, and does not permit an analytic 
solution. It also contains a separation constant, arising in the separation of variables solution of the convected 
wave equation, that acts as the eigenvalue for the unsteady acoustic modes. Numerical efforts to determine 
the modes require a discretization in the radial direction and a subsequent calculation of the eigenvalues and 
eigenvectors of the discretized system. 

Early work in the modal analysis of nonuniform mean flows concentrated on integrating the second-order 
ODE for the radial behavior, generally through numerical means. In the cases of flows with mean axial 
shear, the relevant ODE was first analyzed for two-dimensional flows by Pridmore-Brown [8]. Subsequent 
investigators integrated the ODE using Runge-Kutta methods [9, 10], iterative methods [11], and Galerkin 
based methods [12]. See [13] for a review up to 1975. Since that time, numerical work has been done in [14] 
and [15], and a theoretical basis for shear flows is given in [16, 17]. 

Less attention has been focussed on ducted, compressible swirling flows in turbomachinery, although 
there is a rich literature [18] concerning rotating incompressible fluids without axial flow components. The 
linearized Euler equations for flows with swirl are presented in [19, 20, 21], and cases with solid body swirl 
have been examined in [22, 23] and [24]. Finally, Wundrow in [25] has examined swirling potential flows, i.e., 
flows containing a free vortex swirl. The swirling flow results have been applied to special classes of flows, 
such as solid body and free vortex swirl. This restriction is relaxed in the present analysis. 

While it is intuitively reasonable to assume that the unsteady behavior in time and the circumferential 
direction can be described by exponentials, it is not as clear that this is an appropriate assumption for 
the axial direction. There are thus significant existence, uniqueness, and completeness questions regarding 
the analysis. These questions have been addressed to some degree by Shankar [26, 10] and Swinbanks [27], 
who suggest that completeness may be achievable. Goldstein [16, 17] described the complete solution for 
transversely sheared mean flows, and demonstrated the existence of a continuum of eigenvalues, some of which 
give rise to acoustic disturbances and some of which do not. These two categories of disturbances are discussed 
to some degree below. It is notable that Case [28, 29] specifically warned against the “naive” approach of 
assuming exponential behavior, which could potentially miss the continuum modes. He recommended that 
the linearized equations be solved using Fourier and Laplace transforms. In doing so, the continua of 
eigenvalues arise naturally. Perhaps even more interesting is the fact that, according to Case [28], the same 
continuum was first observed by Rayleigh [30] in 1913. 

In this report, the flow inside a constant radius annular or cylindrical duct with acoustically lined walls is 
examined using the linearized unsteady Euler equations. Any unsteady disturbances are presumed to have 
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exponential dependence in the axial and circumferential directions as well as time. The radial behavior is then 
found by discretizing the equations on a Chebyshev-Gauss-Lobatto grid and analyzing the discrete system 
using a standard eigensystem package. In the nonuniform mean flows studied, the unsteady disturbances 
fail into two groups: a discrete set analogous to the acoustic modes in uniform mean flows and a continuous 
band of “converted” modes. 
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Figure 2.1: Annular duct geometry. 

2. Aerodynamic Model 

The physical system consists of a cylindrical or annular duct (see Fig.2.1) with or without acoustically 
lined walls. It is assumed throughout this section that the mean flow is isentropic, the steady velocity has 
no radial component, and that all mean flow quantities are at most dependent upon radial position. The 
governing equations are the unsteady Euler equations, which can be linearized about a nonlinear mean flow. 

2.1 Mean Flow 

The mean flow (and therefore not time-dependent) entropy, momentum, and continuity equations 

VS = 0 

(V-V)V = -p~ l VP 

and 

v - (Vp) = 0 

where 5, V, P, and p represent the entropy, velocity, pressure, and density in the mean flow. 

The mean flow is taken to be isentropic and axisymmetric, so that it has the general form 

V (r, 0, x) = V x (r)e x + V e (r)e e (2.4) 

where e x and e* are unit vectors in the axial and circumferential directions, respectively. 

In addition to handling general mean flows satisfying (2.4), several special flow cases are examined in the 
Results section below. Chief among these are (a) axial shear flow, (b) solid body swirl, and (c) free vortex 
swirl. Mean flows containing only axial shear, i.e., V = V x {r)e x , have been examined by many authors and 
are thus used for validation of the present analysis. For such mean flows, the mean pressure, density, and, 
therefore, speed of sound, are all constant across the duct! 

If the mean flow contains a swirl component, the mean pressure, mean density, and mean speed of sound 
are no longer uniform. Integration of the radial momentum equation leads to 

P = f/^d f («) 

where f is the radial coordinate normalized by the tip diameter ry. Thus the reference point is the duct 
location, f = 1. The isentropic condition VS = 0 implies that VP = A 2 Vp and is used to determine the 
mean flow density. With swirl, the mean flow speed of sound can be determined from A 2 = kP/ p where k 


are 

( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 
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is the ratio of specific heats. Defining A = A(r)/A(l) = A(r)/A T where A T is the speed of sound at the tip 
implies that 

A(r) = exp ^^-2“ J dfj (2.6) 

For the special cases of solid body and free vortex swirl, it is easier to determine A by taking the derivative 
of A 2 = nP/p and substituting in the expressions for P and p given above. This produces 


dA = ( J1 -i)yl 
dr 2rA 


(2.7) 


whose solution is 

A 2 {r)= A 2 t -(k-\ )J h . dr (2.8) 

which is normalized by At subsequently. For solid body swirl, free vortex swirl, and the two together, the 
tangential velocity is given by 


V*, b 

V«,fv 

^,sb,fv 


and the normalized speed 

of sound A is found to be 

% b(f) 

- >-(V) a, < i -' a > 

Al(f) 

■ -(VMM 

^sb,fv(^) 

■ -O; 1 ) 

n 2 (i-f 2 )-r 2 (i-^-) +2(#s-i)flfiogf 



(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 


where Cl = CItt/At and f = T /{rr A t). 

Note that Crocco’s theorem implies that for an irrotational flow (which requires V x = constant and 
Vj) = r/r) the stagnation enthalpy is a constant, so the ratio A(r x )/A(T- 2 ) can be calculated from the 
isentropic flow relations between tiny two radial locations r x and r 2 . The mean flow is assumed to be 
isentropic, so the mean flow pressure and density can be expressed in terms of A as 


p(f) = A 2 ^ K ~ 1 \ and P(f) = A 4 /** -1 ) 


(2.15) 


where the pressure, density, speed of sound, and radius have all been nondimensionalized by their values at 
the tip radius, ry. 

According to Kerrebrock [20, 21], solid body swirl can be used to represent the flow in a behind a high- 
work blade row, and adding free vortex swirl models the flow behind a rotor. The case of purely free vortex 
swirl is also of particular interest, since it is irrotational and can thus be analyzed using the considerable 
machinery developed for potential flows [31, 25]. 


2.2 Unsteady Flow 

The first order unsteady entropy, momentum, and continuity equations are 


Ds 

Dt ~ ° 


Dx 


— -s*(V V)V + (v • V)V = -V(p _ 1 p) 


and 


(2.16) 

(2.17) 
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^ + A 2 (v-V)p + pA 2 Vv = 0 (2.18) 

where D/Dt is the convective derivative operator following the mean flow and tildes indicate perturbation 
quantities. The present analysis assumes isentropic flow, so s = 0 which implies p = A 2 p. 

Axial shear flows 


In an axial shear flow, the governing equations (2.17,2.18) can be expressed in cylindrical coordinates as 

1 dp 


dv 3 
~dt 

i 

pA 2 


dv r , T/ dv r 

dt x dx 

dv$ . T/ dve 
dt x dx 

+ V x ^ + Vlv r = 


dx 


p dr 
1 dp 
pr d6 
1 dp 
pdx 


?E + v x ^) 

dt X dx) 


dv T v r 1 dve dv x _ n 

dr + r r dO dx 


(2.19) 

(2.20) 
(2.21) 
( 2 . 22 ) 


If one then assumes that the perturbation quantities p,v r >v$ y and v x all have the exponential dependence 

/(r, 8 , x, t) = /(r)c i ^* +mtf “ wt) (2.23) 

for / = p, v r , v $ , v x , then the velocity components can be expressed in terms of the pressure and its derivative 
via 

-i dp 


v r = 


VQ = 


V x = 


(k - M x 7) dr 
m 

f (Jk — M x 7)^ 
-M z dp 




(Jfc - M x 7) 2 dr (k - M x 7) 

where A: = urr/Ar and 7 = fc z r T . When substituted into the continuity equation, these give 


dPp 

dr 2 


+ 


2A^7 


dp 

(k-M x 7)] df 




m 


P = 0 


(2.24) 

(2.25) 

(2.26) 

(2.27) 


as the second order ODE for the pressure. 

Previous attempts to determine the radial behavior of the unsteady quantities in an axial shear flow 
concentrated on solving this ODE numerically. This was examined in the present work, and is discussed in 
Section 3.5 below, it should be mentioned that this approach can miss entire families of perturbations. As 
illustrated in Goldstein [32], the linearized momentum and continuity partial differential equations governing 
the behavior of parallel shear flows can be combined into a single partial differential equation by taking the 
divergence of the momentum equation and subtracting the material derivative of the continuity equation, to 

get 

A-H-R - V 2 p = 25V"— (2.28) 

A 2 Dt 2 P P x dx K 1 

where V x ' indicates the radial derivative of the axial mean flow velocity. A single equation for pressure can 
be found by taking the material derivative of both sides and using the radial momentum equation, which 
gives 

1 D 3 p D /—2~\ , ot ,. 9 2 p 


^ d «.-£< v! « +2 ^ = ° 

The exponential assumption (2.23) applied to this equation gives 

’I 2M^7 


(k - M x 7) 


[r ( k-M x ^)\df 


dp 

-4z + 


(k 


- M x 7) 2 - - 7 2 P j = 0 


(2.29) 


(2.30) 
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Case [28, 29] and Goldstein [16, 17] each point out that equations of the form 

x/(x) = 0 (2.31) 

where f(x) is a differential operator have as solutions x — 0, f(x) = 0, and f(x) = <5(x), where <5 is the 
Kronecker delta. In the case of (2.30), this becomes 

A L\p] = 0 (2.32) 

where A = k — 7 M x and L\p] represents the second order ODE for pressure. There are thus three sets of 

solutions to the shear flow eigenvalue problem: (1) the case A(r) = 0 everywhere across the duct, which gives 

rise to a continuum of disturbances; (2) the case L[p(F)] = 0, which produces a discrete set of eigenvalues 
7fc arising from the ODE; and (3) the case L\p(f)] = <5(A = 0), which also consists of a continuum of 
eigenvalues. The circumferential momentum equation with the exponential assumption reduces to \v$ = 
— im/(pr)p , indicating that the disturbances in set (1) have no pressure associated with them. Set (2) are the 
wavenumbers analogous to the discrete acoustic disturbances in uniform axial flow. The right hand side of 
the condition in set (3) is nonzero only at points where the phase velocity of the disturbance v# — A;/ 7 = M x 
in the duct. This latter set acts as a source for unsteady acoustic disturbances in the flow. 

Flows with mean axial shear and swirl 


The linearized perturbation equations for swirling flows in cylindrical coordinates are 

dir ^ V S dv r . .. dv r 2 Vg . 1 dp V? . 

at r o0 ox r p or prA 2 

die Ve dv$ v dvg fV) . dV e \ . 1 dp 

dt r d6 x dx \ r dr / r pr d6 

dvx_ ,Yi^L a,v — 4.—- - --&E 

dt r d6 1 dx dr Vr pdx 

1 (dp + v + Yl ~ 4- dir ^ | 1 _ 

pA 2 V dt r dO x dx) rA 2 Vr dr r r dO dx 


(2.33) 

(2.34) 

(2.35) 

(2.36) 


where .4(r) is the local speed of sound in the mean flow. 

The same procedure that resulted in a single partial differential equation for axial shear can be applied 
to the swirl flow equations, but the results are not as satisfactory. Subtracting the material derivative of the 
continuity equation from the divergence of the momentum equation gives 


R(±-Ei)- V > 

Dt \pA 2 Dt) 



D_ 

Dt 





ldtv\ 
r dO) 


Uv, 


dig dVg _ \ 

-ar + sr v ‘) 


(2.37) 


Unlike the shear flow case, it is does not appear that going to higher order here can result in a single equation 
for p. In view of the simplicity of the discretized first-order system technique detailed below, the general 
PDE has not been pursued. 

After substituting in the exponential form (2.23) and using M x = V x /A, Me = Ve/A, v r = v r /A , v & = 
ve/A , v x = v x /A, and p = p/(pA 2 ), the linearized equations become 


i(i--Mg-jM x )v r + *M e v s = ^ + (/C . '± Mlp (2.38) 

\A r ) r dr r 

. / k 77i _ _ _ \ Me dM e (#c — 1) _ _ im _ /rt 

-* (l - J Me ~ ^ Mx J ve + [~ + “df" + _ 2^ Me \ Vr = _ T P (2-39) 

-i ^ - yMj - 7 v x + €r = ( 2 - 4 0) 

+ ^+[^Ml + l\v r + fv 0 + ijv x = O ( 2 . 41 ) 
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where the mean flow relationships 


(« - 
2 f 


(2.43) 


have been used. Letting 


A = t--M e -iM x 
A r 

M e dM e k - 1 Mf 
f + df + 2 f 

i/ = -A 2 + \Men 


(2.45) 

(2.46) 


then results in expressions for the perturbation velocity components in terms of the pressure and its radial 
derivative 


v T = \-2i- + (k - l)\M 0 ] 

vr L r -1 


. A dp 
p-hz- — 
v dr 


(2.47) 


m\ k — 1 / A 2 \ 1 f / A 2 \ dp 

v$ = ["7I7 + — ( 1 + v)H P+ 2^l 

[7 (1 - k)M1M x ( M e ( 2m 

“■ - [b ( 1 + — M » Ar J + ^(- 


i/r V Ar 


IF p 


1 \ K — 1 


v 2 r 


M X M$ + 1 ^ 


(2.49) 


The ODE for pressure then becomes 

d?p Jl m( i/ + A 2 ) (1 + iMg 


dr 2 \ f 




2 zml 7 [ (« - dM x ] 1 dp 


iM$ [2 im 


Ar I r 


2 f Af L V 7 r J A L 2 f 

B [*■ - - t 2 ] +¥] + ¥ (* ■ + *Y =&■ 

rrdMg 1 m 2 7 2 7 f(l , (1 - nfM^Mx 


t + (*- DA^r|-^- + y + i 


Mg / ... dM x 2m\ll 


The linearized perturbation equations can thus be written in matrix form as 

Ax = ABx 


4*1 1 dp 

f J dr 


where 


-<3-» 

1 diWg (k - 1) 3 

-Mg + -rb+ — - Mg 
r dr 2 r 


dr 2 r 

d 1 (/t + 1) ,, 2 

-= + t + — - Mg 

dr r 2 r 


-~_Mg 

r 


. ( k m n 

- (3 - 


7 k m m 
-* t - ~ M <> 
V A r 


± + 


■ ( k 171 HJT 

~t 7 “ —Mg 
\ A r 
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and A = —17 


(2.53) 


B = 


' M* 

0 

0 

0 

[ €r 

0 

M x 

0 

0 

1 Vg 

0 

0 

M x 

1 

’ X= ft 

0 

0 

1 

M x 

l P ; 


The matrices A and B are then dependent only upon the radial coordinate. Note that the radial derivative 
operator, d/dr , is retained in the coefficient matrices. 

Matrix B in (2.51) has an analytical inverse provided its determinant 


detB = Ml{M x + 1 )(M X - 1) 


(2.54) 


is nonzero, Thus B is singular whenever M x = 0 or M x = ±1 in the flow. If the mean flow is presumed to 
be subsonic, then the only difficulty that arises is when a boundary layer profile is specified and the mean 
flow Mach number is zero at the wall(s). If B is nonsingular, the generalized eigenvalue problem (2.51) can 
be reduced to a standard eigenvalue problem using the analytical inverse of B. Singularity of B does not 
imply that the generalized eigenvalue problem (2.51) cannot be solved, but it would complicate the solution 
of the eigenvalue problem B _1 Ax = Ax. Consequently, the numerical algorithm is kept in the generalized 
eigensystem form. 

Note also that there are several terms in the A and B matrices that are singular at f — 0, i.e., in the 
case of a hollow cylinder. In the analysis, these terms are evaluated by using l’Hopital’s rule as f 0. 


2.3 Boundary Conditions 

The hard wall boundary conditions applied to the system come from the requirement that there be no 
flow into the walls, i.e., 

V r (0) =V r (l)=0. (2.55) 

where a = ru fr t is the hub-to-tip radius ratio for the duct. 

For acoustically lined ducts, the wall boundary condition used by Unruh and Eversman [12] in sheared 
axial flows is 

Ml) = 77t[i-|m(1)]p(1) (2.56) 

v r (a) = VH [l -lM(a)]p(a) (2.57) 

where tjt and t\h are the complex admit tances of the lin ers along the outer and inner duct walls and 
M(l) = 1) + Mq { 1) and M(a) = y/M£(a) 4- M# (a) are the local Mach numbers at those locations. 

These liner boundary conditions are based on the continuity of particle displacement. It will be presumed 
in this report that the same liner boundary conditions can be used in more general, swirling flows, provided 
the correct local Mach numbers are used. 
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3. Numerical Methods 


The system of equations (2.51) is in the form of a nonhermetian generalized eigenvalue problem. The 
eigenvalues are related to the axial wavenumbers of any unsteady disturbances via A = — 17 , and the eigenvec- 
tors contain the radial behavior of v r , ve, v x , and p. In order to determine the eigenvalues and eigenvectors, 
it is necessary to discretize the equations in the radial direction. Two categories of discretizations have been 
implemented: a spectral technique based on the Chebyshev-Gauss-Lobatto points and a finite difference 
method using both second and fourth order differences. Each case will be presented in turn. 


3.1 Spectral derivative matrix 

For any discretization, a derivative matrix is required, such that 


= S’ 

(3-1) 

for an arbitrary discretized function fj = f(xj). The spectral grid uses 
points [33], which satisfy 

the Chebyshev-Gauss-Lobatto 

7T7 

Xj = cos— , j = 0, ... ,7V 

(3-2) 

and aie the extrema of the A^th order Chebyshev polynomial 


Tpt(x) = cos (N cos -1 x ) . 

(3-3) 


Let f(x ) be a smooth function of x on the interval [- 1 , 1], Then f(x) is interpolated by constructing the 
Nth order interpolation polynomial gj(x) such that gj(xk) = and 


N 


M* ) = 53 h i9i( x ) 
j = 0 

where h(x) is the polynomial of degree N and hj = /(xj), j = 0, . . . , N. It can be shown that 


where 


g . — / 2 j = 0,N 
J_ \l 1 < 3 < N — ! 


(3-4) 

(3.5) 

(3-6) 


The derivative of f(x) at the collocation points is then computed by talking the analytic derivative of gj{x) 
and evaluating it at the collocation points, which produces (D^)kj = g’j{xk) where 


r c,(-iy + i 


( Dn)u = 


CjiXl-Xj) 
~ x J 

2(1 ~ x j) 

2N 2 + 1 

6 

2N 2 + 1 


1*3 




1=3 = 1 

l = j=N 


(3.7) 
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is the required derivative matrix [33]. It is shown in [34] that the most accurate method of formulating this 
matrix is to use 


(Dn)ij = < 


- q (- i) t+i 

2£j sin (fc + j) sin ( k — j) 

-cos(jfrQ 

2sin 2 (^0 
2N 2 + 1 


1*3 

l<l = j<N-l 

1=3 = 1 


for l = 1, . . . , N/ 2, j = 1 ,N and then use 


(3.8) 


{D N )ij = l = ~ + 1 , . . . , N 


(3-9) 


for the rest of the matrix. 

A further improvement on this technique is given in [35]. The idea there is to transform the grid such 
that 

sin -1 (ax) 


£ = s(x,a) = 


sin 1 a 


(3.10) 


where a is a parameter, so that the derivative of the function /(x) is now given by 


# _ d^dx 

d£~ dx<% 


where 


dx 


sin 


-l 


s'(x,a) a 

Then the derivative of /(x) at the collocation points is given by 

f'(x) = MD N f{x) 


— y/l — (ax) 2 . 


where the diagonal matrix M has elements 

M n 


1 

s'(xi,a) 


(3.11) 

(3.12) 


(3.13) 


(3.14) 


The new derivative matrix is then defined to be MDn, which shall be referred to as D\. It is shown in [36] 
that am appropriate choice for a is 

a = sech (3.15) 

where e is roughly machine zero 1 . This has been implemented into the code and yields excellent results. 
With this modification, less than three points per wavelength are required to resolve a given eigenmode. 

The derivative matrix Ds is therefore a full matrix, so any blocks in (2.51) containing d/dr will be full. 
The remaining blocks are diagonal, where the radial distribution of the quantity is distributed along the 
diagonal. The resulting eigenvalue problem is complex, nonhermetian, and 4 N x 4 N where N is the number 
of grid points in the radial direction. 

The Chebyshev-Gauss-Lobatto points (3.2) are shown in Fig. 3.1. The grid points are thus concentrated 
at the edges of the domain, which is appropriate for mean flows with boundary layers at the walls. Note 
that the grid points in (3.2) range from x 0 = lto:rjv = -l. The grid is reversed in the present analysis, so 
that the transformation from each system to the other is given by 


r = 


1 - a 
2 



and x 


2 f 

1 — a 



(3.16) 


1 On the IBM RS6000 that the bulk of the results were run on, double precision machine zero is approximately 1.11 x 10 -16 
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a r k 1 


Figure 3.1: Chebyshev-Gauss-Lobatto grid, and mapping from [—1,1] to [<r, 1]. All angles Adi = A6j, and 
the xj grid points are the cosines of the equally-spaced angles. 

The spectral derivative D expressed in Eq. (3.7) was derived based on Lagrangian interpolation polynomials 
in spectral space, i.e., the range [-1, 1]. In the present numerical approach, all calculations other than the 
spectral derivative are done in physical space, so it is easiest to modify D in (2.51) to be 

D. = ^ D (3.17) 

which is what is actually implemented in the code. 

The solution technique is thus a collocation scheme which requires that the system equations be enforced 
at each radial station. Application of the boundary conditions (2.55) is done by replacing the equations at 
r = a and r = 1 by the discrete forms of the boundary conditions. Soft-walled boundary conditions modeling 
the effects of acoustic liners are introduced using a complex admittance coefficient r). 

3.2 Finite difference derivative matrix 

Second and fourth order finite difference derivative matrices were developed for comparison to the spectral 
solutions. Both assume a uniform grid in the x direction. Then, for the formally second order finite difference 
derivative on a uniform matrix, the D matrix defined in Eq.(3.1) above becomes 

-3 4 -1 

-10 1 


-1 0 1 
1 -4 3 



(3.18) 
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and the formally fourth order derivative matrix 

is given by [38], [39] 





' -25 

48 

-36 

16 

-3 





-3 

-10 

18 

-6 

1 





1 

-8 

0 

8 

-1 




Di ' N = 12A* 





1 -8 

0 

8 

-1 






-1 6 

-18 

10 

3 






3 -16 

36 

-48 

25 


(3.19) 


Performance of the finite difference and spectral derivatives is examined in the Results section below. 


3.3 Eigenmode Calculation 


After applying the boundary conditions, the discretized matrices are analyzed using eigensystem routines 
in the linear algebra library LAPACK 2 [40]. The routines use the QZ algorithm to calculate the eigenvalues 
and eigenvectors for a pair of complex, nonsymmetric matrices A and B such that Ax = ABx. 

The QZ algorithm is an implementation of the commonly known QR algorithm to complex, non-Hermetian 
matrices. In brief, the QR algorithm operates as follows [41]: 

Suppose a matrix C is given, and it is desired to determine its eigenvalues and eigenvectors. The Gram- 
Schmidt process is first used to make the columns of c orthonormal, which involves subtracting off the 
projection of the (A: + l)th column of A from the preceding k columns, and normalizing the result. This 
creates an orthonormal matrix Q and an upper triangular matrix R. 

The QR factorization then takes the factorization of C into matrices Q and R and creates a new matrix 
Ci by multiplying in the opposite order: 

Ci = RQ = Q- X QRQ = Q -1 CQ (3 20) 


This matrix has the same eigenvalues as C, since if Cix = Q _1 CQx = Ax, then CQx = AQx. The 
eigenvalues are the same, but the eigenvectors have been scaled by Q. A new QR factorization of Ci gives Qi 
and Ri, and then C 2 = R 1 Q 1 is formed. This process continues iteratively, preserving eigenvalues at every 
step. The sequence C* approaches a diagonal or upper triangular form, and the diagonal entries approach 
the eigenvalues. Shifts are also employed to help convergence. Also, C is preprocessed into a tridiagonal 
form. Various other steps are taken to adapt the algorithm to the generalized eigenvalue problem, and to 
make it more numerically stable and efficient. 

Equation (2.51) is in the form Ax = ABx. Applying the hard wall boundary condition v r = 0 at the hub 
and the tip radius places zeros on the diagonal of matrix B. Wu, et al. [42] suggest that one approach to 
avoid the mathematical difficulty associated with this matrix singularity is map the complex A plane onto a 
disk via the linear transformation 


A = 


r 4- 1 

T — 1 


(3-21) 


The problem then becomes the calculation of the eigenvalues r of the system 


— (A 4- B)x = r(B — A)x 


(3.22) 


The least stable mode then corresponds to the eigenvalue lying closest to the unit circle. Numerical experi- 
ments were run both using this system of equations as well as the original system (2.51), without significant 
differences. 


3.4 Spectral eigenmode evaluation 

When the eigenvalues are determined using the spectral technique described above, the resulting eigen- 
vectors are determined at the Gauss-Lobatto points. The eigenmodes can be evaluated at points in between 
the Gauss-Lobatto points by employing the appropriate Chebyshev interpolation matrix. 

2 The LAPACK source distribution is available by anonymous ftp from netlib.att.com. 
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The interpolation matrix is based on the interpolation polynomial (3.5) at the Gauss-Lobatto points so 
that the value of the eigenmodes at x ^ Xj (where x 3 are the Gauss-Lobatto points) is given by 


f(x) = 


(i - z 2 ros) 
N 2 




Cj(x-Xj) 


(3.23) 


3.5 ODE integration 


As remarked above, calculations of the nonuniform flow eigenvalues in the literature have focussed pri- 
marily on integrating the second-order, linear, variable-coefficient ordinary differential equation for the radial 
behavior of the pressure. This approach was therefore examined in the present work, though ultimately not 
used. This section documents the reasons for doing so and the lessons learned. 

For the shear flow equation (2.27) requires the solution of a two-point boundary value problem in one 
unknown parameter 


(3-24) 


where A = k - M x 7, with the boundary conditions p'{c) = p'(l) = 0. 

To formulate the problem numerically, the equation was transformed to three first order equations in the 
variables p\ = p, P 2 — p', Pz = 7: 


Pi = 
P2 = 
Ps = 


P 2 




j\ 

~^~ Pa ) 


P 1 


0 


(3.25) 

(3.26) 

(3.27) 


with boundary conditions piicr) — P 2 O) = 0* Note that (1) A contains p 3 in it, (2) the first order system 
is nonlinear, and (3) the axial wavenumber has been designated the eigenvalue, since it is unknown and 
expected to be a constant thus satisfying P3 = 0. For the swirl problem, the eigenvalue also appears in the 
boundary conditions. 

This problem can be solved using standard boundary value techniques, such as shooting algorithms or 
relaxation algorithms. In the former, one boundary is fixed (say f = <j) and the equation is integrated to 
the other boundary r = 1. In general the solution will not satisfy the boundary condition at f = 1, so the 
error is used to correct the approximation at the first boundary and the calculation is done again. In the 
latter, relaxation, technique, an initial guess is made for the solution that satisfies the conditions at both 
boundaries but not necessarily in between. This again provides a residual error, which is used iteratively to 
drive the solution to convergence. 

An effort was applied to implement the relaxation algorithm from [43] to solve the uniform, shear, and 
swirl flow ODEs. As an initial condition, the uniform flow radial acoustic modes, calculated by an alternate 
method, were used. 

The algorithm was found to converge relatively well for the discrete set of modes analogous to the 
uniform flow acoustic modes. Attempts were made to map out the domains of attraction of the eigenvalues, 
with mixed results. One significant problem arose, however, regarding the continuum of eigenvalues. This 
continuum occurs when the parameter A defined above is zero. As it approaches zero in the shear flow 
equation (3.24), the coefficient of dp/ dr containing 1/A grew rapidly. The condition defining the existence of 
an eigenvalue is that the ODE have a nontrivial solution at those values of the parameter. In the vicinity of 
A « 0, the entire equation collapsed to a constant times dp/ dr equals zero, which has the nontrivial solution 
p(f) = constant / 0. Consequently, any time the path of integration taken by the iteration algorithm passed 
near A = 0, the iteration would stop and claim convergence. 

In principle this problem can be avoided by doing an asymptotic expansion of the equation in the 
vicinity of A = 0 and proceding from there. While this was viewed as possible, the idea of applying the 
straightforward spectral technique described above on the original system equations proved to be far simpler, 
and also promised to yield the entire set of eigenvalues in one calculation rather than over several dozen 
runs. The results presented below are thus taken from the spectral technique, though the ODE integration 
was used to validate some of the sheared axial flow results. 
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4. Results 


4.1 Validation 


Validation of the discretized eigensystem analysis was performed by comparing to the analytical solutions 
available [44] for uniform mean axial flows and with published results for flows in lined cylindrical and 
annular ducts containing mean axial shear. For a uniform axial flow, the axial wavenumbers of any acoustic 
disturbances are given by 

Tacoustic = ~ 1 _ ^ 1 - M 2 \A 2 + (^*1) 

where K mti is a separation constant for the mth circumferential and fi th radial mode in the solution of the 
convected wave equation. The constants are determined using the boundary conditions. They are 
countably infinite and can be ordered by magnitued. The cut-off frequency is given by 

A = fc 2 - {Ml - 1)/4 m = 0 -»■ fccut-off = y/l - Ml (4.2) 


which gives for the wavenumber 


'Tcut— off — 


kM x 
AT 2 - 1 


(4.3) 


When A is positive, the associated acoustic disturbances propagate unattenuated in the axial direction 
upstream and downstream from their point of origin. Otherwise the acoustic disturbances attenuate expo- 
nentially. The direction of propagation for the unattenuated waves is determined from the axial component 
of their group velocity, 

dk , 7 _ . .... 

- ■ — 4- M x . (4.4) 


u g,x 


d 7 ± v /t 2 + «^ 


The axial wavenumbers for any entropic and/or vortical disturbances, which are purely convected by the 
mean flow, are given by 

= W, (4 ' 5) 

As part of the validation of the numerical techniques, a version of the present analysis applicable to 
lined, two-dimensional ducts was devised. Results for a two-dimensional duct using a uniform mean flow 
at Af* = —0.5, a reduced frequency k = —6.0, and liners with admittance rj = 0.72 4 0.42i are shown in 
Table 4.1, where they are compared with those from [45]. The results taken from [45] are the results of a 
high-order Runge-Kutta scheme used in that paper. For the present analysis, 32 radial points were spread 
across the domain. Note that the two sets of results are essentially identical. 

Table 4.2 uses data from the same reference, but for a duct with a uniform core flow from y = 0 to 
y = 0.8, and a linear Mach number profile from y = 0.8 to y = 1, where the Mach number at the boundary 
was set to zero. The core Mach number is M x = 0.3, the reduced frequency in this case is k = —5.0, and the 
admittance is 77 = 0.1607 4 0.4463i. Once again, the differences are minimal. 

Finally, comparison results for the flow in a lined, hollow cylinder containing a uniform flow are given 
in Table 4.3, where the comparison results are from [15]. This case has m = 2, k = — 1, M x = 0.5, and 
r) = 0.72 4 0.42i. Once again, the results are in good agreement. 

The data from [10] for a cylindrical duct with a shear profile 


V x (r) = Vb(l - r) 1/7 


(4.6) 


axe shown in Table 4.4 for Mo = 0.3 and k = 20. Table 4.5 gives results for an annulus with a = 0.85714 
and k = 70, with the Mach distribution 


V x (r) = V 0 



1 

( l -*) 2 


w 


+ 1 - 



(4.7) 
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Table 4.1: 2D, Uniform Flow with Liner 


7 * 

[45] 

Present 

7o + 

1.964 -0.0037 

1.9636 - 0.0033? 

7i + 

1.622-0.051t 

1.6218 - 0.05067 

7 2 + 

0.979 - 0.736i 

0.9788 - 0.73567 

it 

0.831 - 1.5007 

0.8307 - 1.49997 

7 4 + 

0.753 - 2.2197 

0.7531 - 2.21917 

it 

0.715 - 2.9017 

0.7150 - 2.9009i 

76 

0.691 - 3.5607 

0.6915 - 3.56007 

it 

0.675 - 4.2067 

0.6752 - 4.20607 

it 

0.663 - 4.844i 

0.6633 - 4.8436i 

it 

0.654 - 5.476i 

0.6543 - 5.4755i 

lo 

-0.655 + 0.0457 

-0.6547 + 0.0448i 

h 

-0.592 + 0.074i 

-0.5921 + 0.0740i 

72 

-0.112 + 0.1527 

-0.1118 + 0.15157 

73 

0.609 + 0.973i 

0.6090 + 0.97337 

74 

0.685 + 1.8207 

0.6851 + 1.8199i 

Is 

0.719 + 2.526i 

0.7192 + 2.5257i 

Is 

0.744 + 3.183i 

0.7437 + 3.1829i 

77 

0.762 + 3.815i 

0.7623 + 3.8154i 

7s 

0.775 + 4.433i 

0.7752 + 4.43307 

7 9 

0.782 + 5.042i 

0.7822 + 5.04207 


Table 4.2: 2D, Shear Flow with Liner 


7^ 

[45] 

Present 

7o" 

0.6428 - 0.02157 

0.6431 - 0.02157 

it 

0.8553 - 0.04157 

0.8547 - 0.04107 

it 

-0.2853 - 0.55587 

-0.2852 - 0.55587 

lo 

-1.4200 + 0.04527 

-1.4196 + 0.04517 

h 

-1.2946 + 0.05967 

-1.2945 + 0.05777 

72 

-0.3995 + 0.56347 

-0.3993 + 0.5633i 


The data in Table 4.6 are for a lined annulus with Mq = 0.3, A; = 30, a = 0.66667, and r] = 0.3 4- O.li. 

The lack of available results for flows with mean swirl has hampered the validation efforts in that area. 
An “internal consistency check” has been done, however, in the following manner. Equations (2.47), (2.48), 
and (2.49) can be used to express the perturbation velocity components in terms of the pressure and its radial 
derivative. The eigenvector x in the system equations (2.51) contains the computed radial distribution of all 
four of the perturbation variables v r , ve,v x , and p. An internal consistency check was done by extracting the 
computed results for pressure from the eigenvectors corresponding to acoustic perturbations and substituting 
them into the analytical relations for the velocity perturbations and comparing to the remaining parts of 
the computed eigenvectors. In all cases examined, the comparisons gave excellent agreement. Since the 
computed pressure is being used, this test confirms that (1) the equations represented in the numerical code 
are the desired equations, and (2) high order agreement is achieved throughout the grid. 

4.2 Accuracy of Numerical Discretization Technique 

To test the relative accuracy of the numerical discretization techniques, a series of uniform flow cases in 
a hard-walled hollow cylinder was run using the parameters M 0 = 0.181684, cr = 0.0, and k = 20.0. This 
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Table 4.3: Cylinder, Uniform Flow with Liner 


it 

[15] 

Present 

To" 

0.620 - 5.014* 

0.6195 - 

5.0139* 

it 

-5.820 - 3.897* 

-5.8195 - 

3.8968* 

it 

0.445 - 9.187* 

0.4453 - 

9.1868* 

it 

0.453 - 13.062* 

0.4539 - 

13.062* 

it 

0.480 - 16.822* 

0.4795 - 

16.822* 

it 

0.503 - 20.531* 

0.5029 - 

20.531* 

it 

0.522 - 24.213* 

0.5220 - 

24.213* 

it 

0.538 - 27.880* 

0.5376 - 

27.880* 

it 

0.550 - 31.537* 

0.5502 - 

31.537* 

it 

0.589 - 49.75* 

0.5891 - 

49.754* 

lo 

0.410+ 1.290* 

0.4101 + 

1.2904* 

h 

1.259 + 6.085* 

1.2595 + 

6.0852* 

72 

1.146 + 9.668* 

1.1457 + 

9.6679* 

7 3 

1.022 + 13.315* 

1.0218 + 

13.315* 

74 

0.943 + 16.977* 

0.9425 + 

16.977* 

7s 

0.891 + 20.635* 

0.8908 + 

20.635* 

76 

0.855 + 24.288* 

0.8549 + 

24.288* 

77 

0.829 + 27.937* 

0.8288 + 

27.937* 

Is 

0.809 + 31.581* 

0.8089 + 

31.581* 

79 

0.755 + 49.77* 

0.7547 + 

49.772* 


Table 4.4: Cylinder, Shear Flow 


it 

Shankar [10] 

16 points 

32 points 

lo 

0.81500 

0.81600 

0.81547 

h 

0.76944 

0.76956 

0.76952 

72 

0.72751 

0.72766 

0.72763 

7s 

0.65329 

0.65342 

0.65346 

7 4 

0.54028 

0.54035 

0.54051 

Is 

0.36933 

0.36937 

0.36975 

7e 

0.06361 

0.06491 

0.06477 

7t 

-0.28313 + 0.48807* 

-0.27836 + 0.48717* 

-0.28277+0.48650* 

78 

-0.28357 + 0.80635* 

-0.28775 + 0.86216* 

-0.28328 + 0.80439* 

79 

-0.28410+1.05658* 

-0.26178+1.0964* 

-0.28368+1.0539* 

7io 

-0.28410+1.27947* 

-0.38211 + 1.6308* 

-0.28403+1.2756* 


corresponded to the case examined in Table 1 of [46]. The output parameter of interest in [46] was K = 7/Jfc, 
which was computed from the second-order differential equation for pressure using a 4th-order accurate 
Runge-Kutta method. 

Figures 4.1, 4.2, 4.3, and 4.4 plot the natural log of the error against the natural log of the number of 
radial points, for the eigenvalue corresponding to the Oth, 1st, 2nd, and 3rd acoustic modes. 

The finite difference approximations use an equally-spaced grid and the derivative matrices given in (3.18) 
and (3.19). The spectral discretization uses the Chebyshev-Gauss-Lobatto points and the derivative matrix 
in (3.7). Each figure also shows a dotted line representing formal second and fourth order accuracy. 

In general, the second and fourth order finite difference approximations follow the formal definitions. 
The spectral accuracy is significantly greater for all numbers of grid points, and is the scheme utilized in the 
presentation of further results below. 
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Table 4.5: Annulus, Shear Flow 


in 

Shankar [10] 

Present 

To" 

0.79293 

0.79353 

7i 

0.75075 

0.75292 

72 

0.57143 

0.57320 

73 

-0.00969 

0.16437 

74 

-0.28733 + 0.732191 

-0.28357 + 0.73425* 

7s 

-0.29118+ 1.21721* 

-0.28622+ 1.2198* 

7e 

-0.29248+ 1.62569* 

-0.28766+1.6281* 

77 

-0.29519 + 2.00221* 

-0.28947 + 2.0055* 

78 

-0.29567 + 2.36511* 

-0.29035 + 2.3683* 

79 

-0.29768 + 2.71665* 

-0.29167 + 2.7209* 

7io 

-0.29776 + 3.06414* 

-0.29243 + 3.0679* 


Table 4.6: Lined Annulus, Shear Flow 

7„ 

Shankar [10] 

Present 

7o 

0.78698 + 0.00400* 

0.78093 + 0.00913* 

7i 

0.73438 + 0.02541* 

0.75079 + 0.03387* 

72 

0.55840 + 0.03148* 

0.57267 + 0.03246* 

7 3 

0.14308 + 0.07638* 

0.16875 + 0.06982* 

74 

-0.23900 + 0.74173* 

-0.23734 + 0.72727* 

7s 

-0.26149+ 1.21973* 

-0.25993 + 1.2120* 

7e 

-0.26996+ 1.62627* 

-0.26860+1.6207* 

7 7 

-0.27669 + 2.00192* 

-0.27468+1.9983* 

7s 

-0.27974 + 2.36439* 

-0.27813 + 2.3612* 

7 9 

-0.28359 + 2.71568* 

-0.28147+2.7139* 

7io 

-0.28502 + 3.06304* 

-0.28361 + 3.0610* 


4.3 Swirling Flows 

The remaining results illustrate the effect of swirl on the axial wavenumbers and radial modes of any 
unsteady disturbances. The baseline case is taken to be a mode with two nodal diameters, m = 2, with axial 
Mach number M x = 0.3, reduced frequency k = 10.0, and hub-to-tip ratio a = 0.25. Swirl is provided by 
specifying the amplitude of a free vortex, i.e., 

Vfl(r) = ~ . (4.8) 

Figure 4.5 shows the real vs the imaginary parts of the axial wavenumbers for the baseline case with no 
swirl. The figure shows the classical result, in that the line to the left of the imaginary axis represents 
the analytic cut-off condition, two families of propagating acoustic modes are distributed symmetrically 
about the cut-off wavenumber, and a doubly-infinite set of attenuated modes exist with real parts equal to 
the cut-off wavenumber, Re{7 cu t-off} = -1-5707. Eigenvalues that are purely real indicate unattenuated 
propagation, while for the complex eigenvalues the imaginary part gives the growth or decay rate, c.f. (2.23) 
and (4.1). Physical arguments are used to eliminate growing acoustic modes in either direction. In addition, 
the eigenvalues corresponding to the convected disturbances all fall on top of one another at the convected 
wavenumber 7 CV = 33.3. 

Free vortex swirl with magnitude T = 0.2 is then added to the mean flow, giving the Mach number 
distribution shown in Figure 4.6. The axial wavenumbers from both the uniform flow case and the swirling 
flow case are shown in Figure 4.7. The effect of the swirl on the discrete sets of modes corresponding to the 
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acoustic modes in uniform mean flow is largest near cut off and tends to vanish for large mode order. One pair 
of the uniform flow acoustic modes has become cut off with swirl present. The isolated uniform flow convected 
mode, however, has split into a set of points covering the continuum range from 7 C v,min = kfM x {\) = 12.0 
to 7c v, max = k/M x {a ) = 32.0, since the nonuniform mean flow implies different velocities at different grid 
points. Adding additional grid points to the calculation merely covers the continuum more densely without 
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Order of Accuracy 



Figure 4.3: Order of accuracy for mode 2. 

Order of Accuracy 



Figure 4.4: Order of accuracy for mode 3. 


changing its minimum or maximum values. 

Figures 4.8 and 4.9 illustrate the real and imaginary distributions of pressure for the «2,o an d « 2,2 modes, 
respectively. The radial acoustic modes are indexed by the number of times their real parts cross zero. It 
is interesting to observe that the unattenuated mode /c 2 ,o is almost purely real, while the exponentially 
decaying mode k; 2 ,2 is complex. 
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Axial Wavenumbers 



Figure 4.5: Uniform mean flow wavenumbers. The discrete set symmetric about the cut-off line correspond 
to acoustic modes propagating upstream and downstream, while modes on the cut-off line are attenuated as 
they move axially. The point labelled “pure convection” is the repeated wave number for all the convected 
modes. 


Mach Distribution 



Figure 4.6: Mach distribution for swirl case. 


The eigensystem analysis computes the complete eigenvector given by x = {v r ,ve ,\ and thus com- 
putes the radial distributions of the radial, circumferential, and axial velocities along with the pressure. 
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Axial Wavenumbers 



10 15 

Re{gamma} 


35 


Figure 4.7: Effect of free vortex swirl. 



Figure 4.8: Cut on pressure mode « 2 ,o for (a) uniform flow and (b) free vortex swirl with T = 0.2. 


Simple algebra can be applied to (2.51) to express the velocity components in terms of pressure and its 
radial derivative, and this provides a convenient check of the consistency of the algorithm. Figures 4.10 and 
show the real and imaginary parts of all the disturbance quantities for modes ^ 2,0 and ^ 2 , 2 , respectively, 
both with and without swirl. For the uniform axial flow cases, the velocity modes do correspond to the 
analytic formulas applied to the acoustic part of the overall eigenmode. The qualitative effect of adding the 
swirl can be seen in Figure 4.11. Note in particular that the slope dp/dr of the acoustic eigenmodes at the 
hub and tip boundaries is not zero in the swirl cases as it is in uniform flow. This situation is particularly 
evident at the hub, where the swirl Mach number is greatest. 


23 






















5. Unresolved Issues 


This section documents issues that have not yet been resolved in the present analysis. 


5.1 Hydrodynamic stability 


A necessary but not sufficient condition for instability for axisymmetric (i.e., m = 0) disturbances in 
the inviscid limit is the vanishing of — : (7^) somewhere in the flow 1 . This condition is analogous to the 
existence of an inflection point in a two-dimensional shear flow profile. Also, if an instability is to exist (also 
for m = 0) then the phase velocity of the disturbance v<p = k/j must be equal to the mean flow velocity at 
some point in the flow. 

In an inviscid fluid in a cylindrical or annular duct, Kelvin’s theorem states that the circulation along 
any ring at radius r will be a constant. For a swirling mean flow, the circulation is given by T — 2tt rV$ — 
constant. The following stability argument is due to Von Karman [49]. Define the quantity 


k(r) = rV$(r) (5.1) 

then the centripetal acceleration at any radius is given by Vgjr = k 2 /r 3 . Consider a ring at radius 77 , having 
velocity V$ % \ = V$(ri) and fci = r{Ve y \, and displace a fluid particle in this ring to a radius r 2 > tv 
T he centripetal force on the displaced particle is then given by 


while the pressure gradient at r 2 is 


/centrip — P 3 


x — 

/pressure — P 3 


(5.2) 


(5.3) 


so the fluid is stable (i.e., the restoring pressure force is greater than the displacing centripetal force) if 

k\ 


4 k\ 

r 3 r 3 

r 2 '2 


(5.4) 


implying that k 2 (r) increases outward. 

For a solid body swirl, V$ = Hr so k = fir 2 . This increases outward for all r, and is thus stable. For a 
free vortex swirl, Ve = T/r so k = I\ This is a constant, which would seem to imply neutral stability. In 
order to have a test case that would give instability, a swirl flow of the form V$ = 0jr 2 was implemented. 
This gives k = /3/r, which decreases outward. 

In order to examine the stability of unsteady perturbations in the flow, a test due to Briggs [50] (and 
used in Tam and Hu [51] and Hu [52]) is used. It involves deformations of contours in the complex plane 
as part of an inversion of a Laplace transform. Letting u be complex, the eigenvalue problem is solved for 
a series of frequencies u = cjr + to;/, where ur is the desired frequency and uj is steadily reduced to zero. 
Any zeros originating in the upper half of the 7 plane are downstream propagating waves, and zeros from 
the lower half of the 7 plane are upstream propagating waves. Any zeros that move across the real 7 axis are 
instability waves. All the others either propagate unattenuated (neutrally stable) or decay exponentially. 

Figures 5.1 and 5.2 illustrate this process applied to the special swirl cases V$ = T/r and V$ = /3/r 2 . The 
axial flow was uniform at M x = 0.2, along with u — 0.5, 0 = 0.2, and m = 2. For both figures, the discrete 
eigenvalues move from their values for k = 10 + 2i to their values at k = 10 smoothly, without crossing the 
real axis. The range of convected modes between Re{7} = 32 to Re{7} = 48 for the free vortex swirl case 
also do not cross the real axis. 

In the contrived case V$ = /3/r 2 , however, the convected modes have split into an oval pattern, which 
(when k is purely real) is symmetric about the real axis. Consequently, when k is given a finite imaginary 
part, part of the oval crosses the read axis indicating instabilities. 

The oval nature of the convected roots has not been verified by previous investigations. A grid resolution 
study appears in Figure 5.3 (note the expanded ordinate scale). While the discrete roots have clearly 
converged, the oval structure is still shrinking toward the real axis as more and more points are added. While 
it is possible that this behavior is an illustration of so-called pseudo-eigenvalues 2 and therefore physically 


1 This condition, attributed to Rayleigh, is discussed in [47] and [48]. 

2 If eigenvalues are defined as A such that ||A - AI|| =0, then the pseudo-eigenvalue A e satisfies 1 1 A — A e I|| < e. 
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Figure 5.2: Stability of V$ = /3/r 2 . 

real, further checks are required. 
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Figure 5.3: Grid independence study for oval structure. 


5.2 Completeness 

The form of the unsteady perturbations assumed in (2.23) is equivalent to taking Fourier transforms in 
the axial and circumferential directions and in time. While considerable evidence in favor of this assumption 
exists for the independent variables 6 and t, the situation is not as clear in the x direction. In fact, recent 
work by Goldstein and Wundrow [25] suggests that the vortical part of the unsteady velocity in a swirling 
potential flow (i.e., free vortex) has an algebraic dependence on x. Tan[53, 54, 55], however, has suggested 
that, the dependence is linear in x near the source of the disturbance but proportional to l/x later. Resolution 
of this issue is critical for completeness, since if the modes have an algebraic dependence on x, the axial 
behavior is neither bounded nor square integrable so the exponential assumption for the axial behavior is 
not valid. 

It is interesting to note that the exponential assumption arises naturally from a separation of variables 
type analysis. The matrix equations are given by (2.51). If the form for the unsteady perturbations is 
generalized to 


tv(r, d, x, t ) = v r {r)X{x)e i( - m6 -^ (5.5) 

v 9 (r,e,x,t) = v e (r)X(x)e i{m6 - ut) (5.6) 

v x (r,6,x,t) = v x (r)X(x)e« me - ut) (5.7) 

p(r, 9,x,t) = p(r)X (x)e^ me -^ (5.8) 


then a matrix differential equation can be written for the axial behavior, X (x), 


where 


Au 



V x 0 0 0 

0 V x 0 0 

0 0 V x 1 

0 0 1 v x 


(5.9) 


(5.10) 
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said 


B = 


i{?V e - w) -2^. 0 

^ + i(*V,- U ) 0 

, d -t 0 i(*V t - U ) 

^ + \ + D i? i(fV t - u ) 



(5.11) 


and u = {v r ,ve,v x ,p} T . This reduces to the form given in (2.51) when -X r /X = A = constant. If (5.9) is 
viewed as an ordinary differential equation in x , then it can be solved in terms of matrix exponentials as 


X(x) = X 0 e- B " lAx 


(5.12) 


This process fails if B is singular, which will occur at no-slip boundaries where M x = 0: 

It should be noted that if the axial wavenumbers in Eqn.(2.23) are assumed a priori to have a radial 
dependence, i.e., 


tv(r, x) = v T {r)e %1 ^ z 


(5.13) 


and likewise for the other unsteady variables, then the radial derivatives 
additional term: 


8v r 

dr 


( 


dv r . d'y 


X ) 


e 


iyx 


appearing in (2.51) produce an 

(5.14) 


This does introduce an algebraic dependence on x, and must be examined further. 

An alternative approach, which takes advantage of the existing theorems on completeness, is to use the 
symmetric form of the nonlinear Euler equations [56]. Subsequent linearization will maintain the symmetric 
form, and the theory states that any normal operator can be both diagonalized and symmetrized, thus giving 
a complete set of eigenvalues. This approach has not yet been pursued, but appears to be viable. 
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6. Conclusions 


The eigenvalues and eigenvectors for nonuniform mean flows in annular and cylindrical ducts have been 
calculated using a discretized form of the first order linearized Euler equations. Three numerical discretization 
techniques have been applied to the first order system; a Chebyshev collocation technique, second-order finite 
differences, and fourth-order finite differences. The resulting complex, nonhermetian, generalized eigenvalue 
problems were solved using the QR algorithm from LAPACK. The analysis has been validated for uniform 
flows and for shear flows with and without acoustic liners. Results have also been presented for annular 
ducts containing free vortex swirl. 

The axial wavenumbers of the nonuniform mean flows are seen to fall into two distinct classes: a discrete 
set corresponding to the acoustic modes in uniform mean flows, and a continuous set along the local convected 
speeds of the mean flow. 

Although the code has been validated for lined annular and cylindrical ducts containing axial shear 
flows, the lack of published results for swirling flows has hampered efforts in this area. Instead, the internal 
consistency of the code has been verified. The spectral discretization technique has been compared to the 
finite difference approximations and found to be significantly more accurate for a given number of radial grid 
points. 

Questions still remain regarding the completeness of the calculated eigensystem, particularly in reference 
to any algebraic dependence in the axial direction of the unsteady disturbances. Efforts to resolve this issue 
will be made, partly through the use of assumed radial dependence of the axial wavenumbers. In addition, 
the stability test advocated by Briggs has been used to examine ostensibly unstable flows, but, other than 
the appearance of an oval shaped region around the convected wavenumbers, no instabilities have been 
detected. Work must also be done to determine whether the oval region itself is an indication of physical 
pseudo-eigenvalues in the system or simply a sensitive detector of subtle errors in the analysis. 
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8. Appendix: The Swirl Code 


The algorithms described above for analyzing axial flows with mean shear and swirl in lined and unlined 
ducts have been implemented into a code named SWIRL. 

The code is written in Fortran 77 with some relatively standard extensions, all of which eventually 
became part of the Fortran 90 standard. Double precision complex data types were used to try to achieve 
the maximum possible numerical accuracy. To facilitate interaction with the user, the “NAMELIST” has 
been employed, in a form compatible with the Fortran 90 standard. It should be noted, however, that the 
specific format of NAMELIST interaction statements varies from system to system. For example, on the 
IBM RS6000 workstation where the code was developed, the user types “ &inputs omega = 10.0 &end” to 
communicate with the program. 

To faciliate both portability and reliability of the program, standard freely-available numerical libraries 
were used wherever possible in the algorithm. These include the numerical linear algebra library LAPACK 
(which also makes use of the Basic Linear Algebra Subroutine, or BLAS, library), and the basic curve fitting 
routines contained within the FITPACK library. Rather than expect the user to download and install these 
libraries, the revelent routines were extracted from them and are included in different subdirectories of the 
SWIRL code distribution. 

8.1 Input Quantities 

A sample input file, input. data, is included in the distribution. The sample file is reproduced here for 
reference: 


& inputs 


mm 

= 2, 

npts = 

32, 

sig 

= 0.50, 

akre = 

10., 

akim = 0 . , 

ir 

= 1, 

rmax = 

0.30, 

slope 

= 0.00, 




is 

= 1. 

angom = 

0.50, 

gam 

* 0.00, 

irepeat 

= 0, 

ifd = 0, 

itest 

= o. 

etahr = 

0.00, 

etahi 

= 0.00, 

etadr = 

0.00, 

etadi = 0.00, 

ed2 

= 0. 

0, ed4 = 

0.0 







&end 

where the various quantities are described in Table 8.1. The sample file thus represents a search for the axial 
wavenumbers for 2-lobed modes of frequency lj = 10 in a hard- walled duct with radius ratio 0.5, containing 
a uniform axial mean flow at mach number M x = 0.3 and a solid body swirl flow with Q — 0.5. Thirty- two 
points are distributed on a Gauss-Lobatto grid across the duct, implying that the resulting wavenumbers 
associated with modes with up to 9 zero crossings are expected to be reproducable. 

The code itself contains a significant number of comments intended as documentation. Of particular 
interest are the subroutines determining the form of the mean flow in the duct, namely, rmach.f, smach.f , 
and sndspd.f. The first, rmach.f , contains options for choosing the form of the axial shear flow. The second, 
smach.f describes options for the mean swirl distribution, and the last, sndspd.f determines the radial 
dependence of the mean speed of sound. Note that a specified distribution can in principle be entered for 
all three, but that the resulting speed of sound distribution is not likely to satisfy the isentropic conditions 
assumed in the analysis unless this was taken into account ahead of time. 

8.2 Output Files and Testing 

Two test output files are included in the distribution: gammas.test and gam.non.test The former contains 
the real and imaginary parts of all the axial wavenumbers calculated by the code. The latter contains only 
those wavenumbers associated with the nonconvected modes (generally acoustic), sorted by number of zero 
crossings in the mode. The code generates the analogous files gammas.dat and gam.nonconv , which should 
be compared to the above, and the file cv.waves.dat which holds the data for the convected waves. The 
wavenumber comparison is most easily carried out using a plotting package like the freely-available program 
GNUPLOT. To facilitate comparisons using this program, any lines in the output files containing nonnumeric 
data begin with an initial “#”. 
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Table 8.1: Description of Input Quantities 


mm 

npts j 
sig 
akre 
akim 

Circumferential mode number 

Number of radial mesh points 

Hub-to-duct radius ratio, a 

Real part of frequency, k = wr^/Ar 

Imag part of frequency (zero except for stability calcs) 

ir 

Axial Mach distribution number (for details see rmach.f) 

0 = Uniform 

1 = Linear shear 

2 — Read from file 

3 = Uniform + sine wave bndry layers of thickness 6 

4 = Uniform + linear bndry layers 

5 = Uniform + 1 /7th power law bndry layers 

6 = Hyperbolic secant profile 

7 = Laminar mean flow 

8 = Wavy sinusoid (for stability calcs) 

9 = Hagen-Poiseuille flow 

rmax 

Maximum axial Mach number 

slope 

slope of linear Mach distribution; 

also used as boundary layer thickness 

is 

Swirl Mach distribution number (for details see smach.f) 

0 = No swirl 

1 = Solid body swirl 

2 = Free vortex swirl 

3 = Solid body and free vortex together 

4 = V$ = 1/r 2 (for stability calcs) 

5 = Read from file 

6 = Constant swirl across the duct 

7 = Trailing line vortex 

angom 

Magnitude of solid body swirl V$ = fir; 

also used as constant for constant swirl 

gam 

Magnitude of free vortex swirl V$ = T/r 

ifd 

Use finite differences for derivatives 

1 = Second order 

2 = Fourth order 

itest 

Perform consistency test on selected modes 

etahr 

Real part of hub liner admittance 

etahi 

Imag part of hub liner admittance 

etadr 

Real part of duct liner admittance 

etadi 

Imag part of duct liner admittance 

ed2 

Second order smoothing for derivatives 

ed4 

Fourth order smoothing for derivatives 
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flows. The physics of swirling flows is examined with reference to hydrodynamic stability and completeness of the 
eigensystem expansions. The effect of assuming exponential dependence in the axial direction is discussed. 
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